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Abstract 

In this work, we develop a novel Galerkin-Ll-POD scheme for the subdifFusion model with a 
Caputo fractional derivative of order a € (0,1) in time, which is often used to describe anomalous 
diffusion processes in heterogeneous media. The nonlocality of the fractional derivative requires 
storing all the solutions from time zero. The proposed scheme is based on continuous piecewise linear 
finite elements, LI time stepping, and proper orthogonal decomposition (POD). By constructing an 
effective reduced-order scheme using problem-adapted basis functions, it can significantly reduce the 
computational complexity and storage requirement. We shall provide a complete error analysis of 
the scheme under realistic regularity assumptions by means of a novel energy argument. Extensive 
numerical experiments are presented to verify the convergence analysis and the efficiency of the 
proposed scheme. 
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1 Introduction 

In this work, we consider the following model initial-boundary value problem for u{x,t): 

dfU — Au = f, in n T >t > Q, 

u = 0, on dfl T > t > 0, (1.1) 

u(0) = V, in n, 

where 17 is a bounded convex polygonal domain in (d = 1, 2,3) with a boundary 917 and n is a given 
function defined on the domain 17 and T > 0 is a fixed value. Here 9f u (0 < a < 1) denotes the left-sided 

Caputo fractional derivative of order a with respect to t and it is defined by (see, e.g. [TSl pp. 91]) 

where r(-) is Euler’s Gamma function defined by r(a:) = s^~^e~‘^ds for a; > 0. 

In recent years, the model (DU) has received much interest in physical modeling, mathematical analysis 
and numerical simulation. The main engine that has fueled these developments is its extraordinary 
capability for describing anomalously slow diffusion processes, in which the mean square variance of 
particle displacements grows sublinearly with time, instead of linear growth for a Gaussian process. 
At a microscopic level, the particle motion is more adequately described by continuous time random 
walk, whose macroscopic counterpart is a differential equation with a fractional derivative in time |24j . 
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Nowadays the model has been successfully employed in many applications, e.g., thermal diffusion in 
fractal domains [26j . ion transport in column experiments [5], and non-Fickian transport in geological 
formation [5], to name just a few. 

Numerically, the presence of the fractional derivative d^u has two important consequences. First, 
the nonlocality in time incurs huge storage requirement as well as much increased computational efforts 
along the evolution of the time. Second, the solution operator has only very limited smoothing property: 
the problem has at best order two smoothing in space EH, and the first derivative in time is usually 
unbounded, cf. Theorem ED in the appendix. These represent the main technical challenges in the de¬ 
velopment and analysis of robust numerical schemes for reliably simulating subdiffusion. The challenges 
are especially severe for “multi-query” applications, e.g., inverse problems and optimal control, where 
repeated solutions of “analogous” forward problems are required, e.g., due to variation in problem pa¬ 
rameters or inputs. To reduce the storage requirement, a number of useful strategies have been proposed, 
e.g., short-memory principle and panel clustering pSlHl ETlf^H] . 

In this work, we shall develop an efficient strategy, called the Galerkin-Ll-POD scheme, for reliably 
simulating the subdiffusion model ini» by coupling the Galerkin finite element method (FEM) with 
proper orthogonal decomposition (POD) to reduce the computational complexity of repeatedly simulating 
subdiffusion, which is important for solving related inverse problems and optimal control. POD is a 
popular model reduction technique, and it has achieved great success in reducing the complexity of 
mathematical models governed by differential equations; see [niiiiMiiisiiiiiii] for a rather incomplete 
list. It is especially attractive in optimal control [niaiiaEn] and parameter inversion uni US]. To the 
best of our knowledge, this work represents the first application of the POD for the subdiffusion model 


(l.I) with a complete error analysis. 


Next we describe the proposed scheme. Let be a shape regular quasi-uniform partition of the 
domain fl, and Xh be the associated continuous piecewise linear finite element space. Meanwhile, we 
discretize the Caputo fractional derivative d^u{t) by the LI approximation d^u{tn) (with a time step 
size t ) [201ES] 


d^u{tn) = 


n—1 


j=0 


r“r(2 — a) 


where the weights {bj} are defined by ( |2.4| ). With the Galerkin FEM in space and LI approximation in 
time, we arrive at the following fully discrete scheme: find G Xh for n = 1,2,... ,N 

+ (VG^", Vv.) = ifitr.),^) G Xh, 

with G Xh being an approximation to the initial data v, where (•,•) denotes the inner prod¬ 
uct. The term involves all solutions preceding the current time step n, indicating the 

computational challenge. In this work, we shall adopt the POD methodology to overcome the challenge. 
Specifically, we take the fully discrete solutions {UJ^}n=o and fractional difference quotients 
as snapshots to generate an optimal orthonormal basis {'4’j}j^i- Since these snapshots are sampled from 
the solution manifold, the POD basis is automatically adapted to the characteristics of the manifold and 
is expected to have good approximation property. Then we employ a Galerkin framework using the POD 
space X™, m < r, spanned by the first m POD basis functions, i.e., find G” G XJ^, n = 1,2..., N such 
that 

(a“t/;^, ip) + {xui, vp) = p) yp g 

with G X™ being an approximation to U°. In the reduced order formulation, the degree of freedom 
is TO, the number of POD basis functions, which is usually much smaller than that of the full Galerkin 
formulation. Hence, it yields an enormous reduction in computational complexity and storage require¬ 
ment. We shall provide a com plete a priori convergence analysis of the scheme. Our main theoretical 

For example for the POD approximation {C/^}))Li generated using the 


3.6 


result is given in Theorem 
Hq{11.)-POT) basis, the following error estimate holds (with ^h = \ log/i|) 


1 

N 


N 


J2\Ht„)-u: 




4«4 






A,), 
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2.3 


for 


where are the descendingly ordered eigenvalues of the correlation matrix K (see Section 

details) under suitable verifiable regularity conditions on the source term / and the initial data v. 

This error estimate consists of three components: spatial error temporal error 0(t“) and 

POD error (X]j=m+iWhile nearly optimal error estimates due to the spatially semidiscrete 
Galerkin FEM is available [131 [El ITT] , it is not the case for temporal discretization by the LI time 
stepping. The LI scheme was first analyzed in [20l [35], where the local truncation error was shown 
to be for twice continuously differentiable (in time) solutions, which is fairly restrictive, cf. 

Remark I A. 1[ Recently some error bounds that are expressed directly in terms of data regularity for the 
homogeneous problem were shown using a generating function approach |14j . however, the analysis does 
not extend straightforwardly to the inhomogeneous case. 

In this work we shall develop a novel energy argument for the LI time stepping to overcome the 
technical challenge in the convergence analysis, which represents the main technical novelty. We shall 
derive optimal error estimates under realistic regularity conditions, and the analysis covers both smooth 
and nonsmooth problem data, cf. Theorem |3.5| Further, the stability result plays an essential role in 
deriving error estimates due to the POD approximation. All the theoretical results are fully confirmed 
by extensive numerical experiments. 

The rest of the paper is organized as follows. In Section]^ we develop an efficient Galerkin-LI-POD 
scheme, and in Section provide a complete error analysis of the scheme. In Section extensive 
numerical experiments for one- and two-dimensional examples are presented to verify the convergence 
analysis. Finally, in an appendix, we briefly discuss the temporal regularity results for problem (l.I). 
Throughout, the notation c, with or without a subscript, denotes a generic constant, which may differ at 
different occurrences, but it is always independent of the solution u, the mesh size h, time step size t, 
and the number m of POD basis functions. 


2 An efficient Galerkin-L 1-POD scheme 

In this section, we develop an efficient numerical scheme, termed as the Galerkin-LI-POD scheme, for 
problem It is based on the following three components: standard Galerkin method with continuous 

piecewise linear finite elements in space, LI approximation in time and proper orthogonal decomposition 
in the snapshot space, which we shall describe separately in the following three subsections. 


2.1 Space discretization by the Galerkin FEM 

First we describe the spatial discretization based on the Galerkin FEM. Let TT be a shape regular 
and quasi-uniform triangulation of the domain D into d-simplexes, known as finite elements and denoted 
by T. Then over the triangulation Tu we define a continuous piecewise linear finite element space by 


Xh = {vh & : n/ilr is a linear function, VT S TT} • 

On the space AT/j, we define the L^(D)-orthogonal projection Ph : —>■ Xh by [Ph’^^x) = {v^x) for 

all X S Xh- Then the semidiscrete Galerkin scheme for problem (I.I) reads: find Uhit) G Xh such that 

(d,^Uh,X) + (Vuh,Vx) = if,x) Vx e t > 0 , (2.1) 


with UhiO) = Vh G Xh- Upon introducing the discrete Laplacian Ah : Xh -G Xh defined by —iAh(p,x) = 
(V:p, Vx) for all ip, x ^ Xh, the semidiscrete scheme (2.1) can be rewritten into 


d^Uhit) + AhUhit) = fhit) t > 0, 


( 2 . 2 ) 


with UhiO) =VhG Xh, fh = Phf and Ah = -Ah- 
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2.2 Time discretization by LI scheme 

For the time discretization, we divide the interval [0, T] into N equally spaced subintervals with a 
time step size t = T/N, and = nr, n = 0,... ,N. Then the LI scheme [5D1 [3S] approximates the 
Caputo fractional derivative d^u(x,tn) by 


d^u{x,tn) = 


r(l -a) — jt. 

’ j=0 b 




ds 


1 


r(l-a) 


E 

i=o 


u{x^ ij+i) — u{x, tj) ^b+i 


(tn - s) °‘ds 




i=o 


T“r(2 — a) 


(2.3) 


where the weights {bj} are given by 

bj = U + iy~°‘ - j = 0,l,...,n-l. (2.4) 

Then the fully discrete scheme reads: given Ujl = Vh & Xh and = Phf{tn) £ Xh^ with Cq = 
r(2 — a), find UJ^ € X^ for n = 1, 2,..., TV such that 


n—1 


{bol + c^T‘^Ah)Uj: = + ^(6,-1 - bj)ur^ + c„t“F, 


h ’ 


(2.5) 




The computational challenge of the fully discrete scheme (2.5) is obvious: To compute the numerical 


solution UJl at t„, the solutions at all preceding time instances are required, as a result of the 

nonlocality of the Caputo fractional derivative d^u. Hence, the computational complexity and storage 
requirement grow linearly as the number n of time steps increases, which poses a significant challenge 
especially for high-dimensional problems and multi-query applications. This naturally motivates the 
development of cheap reduced order models by the POD methodology so as to reduce the effective degree 
of freedom. 


2.3 Galerkin-Ll-POD scheme 

Now we develop an efficient Galerkin approximation scheme based on proper orthogonal decompo¬ 
sition (POD) to circumvent the challenge. We shall first describe the general framework of the POD 
methodology, and then discuss its application to the subdiffusion equation. 

POD is a powerful model reduction technique for complex models, especially time/parameter depen¬ 
dent partial differential equations. It resides on the empirical observation that despite the large apparent 
dimensionality of the solution space (e.g., the degree of freedom of the finite element approximation), 
the solution actually lives on an effectively much lower dimensional (possibly highly nonlinear) manifold. 
POD constructs a problem adapted basis for efficiently approximating the manifold using samples from 
the manifold, often known as “snapshots”, which can be either solutions at different time instances, differ¬ 
ent parameter values, or samples generated using relevant physical experiments. The POD basis functions 
are then employed within either a Galerkin or Petrov-Galerkin framework to generate a reduced-order 
model. 

Now we recall the general framework of POD. Let X be a real Hilbert space endowed with an inner 
product (b ■)x and norm || • ||x- Now for iV € N, let {yn\n=i C X be an ensemble of snapshots and at 
least one of them is assumed to be nonzero. Then we set if = spanjyi, j/ 2 , Vn} C X. Let dim(il) = r 
and let be an orthonormal basis of the snapshot space it. Then any element 2 /„ can be written 

as 

r 

yn^'^{yn,'tpj)x'fpj, n = l,2,...,N. 

T=i 
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POD chooses an orthonormal basis {iljj}T=i for 1 < m < r to minimize the following ensemble average: 


1 


N 


VElls'" 


'j)xi’j\\x- 


( 2 . 6 ) 


A solution of problem (2.6) is called a POD-basis of rank m. This optimization problem is related to the 
jTVxAf corresponding to the snapshots {yn}n=iy which is defined by 


correlation matrix K € 


= N- 




X, 


i,j = 


(2.7) 


By its very construction, the matrix K is symmetric positive semidefinite, and its eigenvectors can be 
chosen to be orthonormal (in the inner product (•r)x)- Further, the number of positive eigenvalues is 
equal to r, the dimensionality of the space it spanned by the snapshots (or equivalently the rank of K). 
The following lemma gives the formula of the POD-basis and the corresponding approximation error 
within the ensemble [31]. 

Lemma 2.1. Let Ai > A 2 > ... > Ar > 0 6e the positive eigenvalues of the eorrelation matrix K and 
ui,..., Ur G be the corresponding orthonormal eigenvectors. Then a POD basis of rank m < r is given 
by 


N 


^0 = 

V n=l 


TlJn 


where (vj)n denotes the n-th component of the eigenvector vj. Moreover, the error is given by 


^'^\\yn-'^{yn,iljj)x'(fj\\x = 


N 


N 






Following the abstract framework, for the subdiffusion model (|l.l|), we choose 2N + 1 snapshots as 

yn = ur\ 

and the fractional difference quotients (FDQs) 


n — 1; 2,..., -f 1, 


n = iV + 2,...,2iV + l. 

The inclusion of FDQs {dlfUf} into the snapshots it is to improve the error estimate below: it allows 
directly bounding the error due to the POD approximation to the fractional derivative term difUf, cf. 
Lemma [2T| In the absence of these FDQs in the snapshots, the error estimate due to POD approximation 

*; see Remark 3.3 for details. The use of difference quotients was 


would involve an additional factor 


first proposed by Kunisch and Volkwein for the standard parabolic equation, and we refer interested 
readers to the recent work |3| for extensive discussions. In this work, we shall follow the work and 
employ the FDQs dlfUJf in the construction of the POD basis. 

In practice, there are several possible choices of the Hilbert space X, and we shall consider two 
popular ones in this work. Our first choice for the POD space is A = Hq(JI) with the inner product 
{u,v)x = (Vu, Vu) for all u, v G i7g(f2). Then the correlation matrix K is given by 


A,,, = {2N + l)-\Vyj,Xy,). 


We denote the corresponding POD basis (called Hq{LI) POD basis) by and the subspace spanned 

by the first m i7g(il)-POD basis functions by XJf, m <r. Then Lemma 
estimate for the POD space Xff 


2.1 


( 2 . 8 ) 
nned 

yields the following error 


1 

2A + 1 


N 


N 




n—0 


i=i 


1=1 


= E 

j=m+l 

(2.9) 
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where are the descendingly ordered eigenvalues of the correlation matrix K. The second choice 

is X = L^(n) with the standard inner product. The correlation matrix K is given by 

k,, = {2N+l)-\y,,yi). ( 2 . 10 ) 

Likewise, we denote the corresponding POD basis (called L^(0)-P0D basis) by and by slightly 

abusing the notation, the subspace spanned by the first m L^(0) POD basis functions by X™. Then in 
view of Lemma 2.1 the POD space X™ satisfies the following error estimate 


1 


2X+ 1 


N 


N 




n—0 




n—l 






where are the descendingly order eigenvalues of the correlation matrix K. 

Next we define the Ritz projection operator ii™ : —>• X^ 


>-h ■ by 


(VR)rx,V(^) = (Vx,v^) Vv^exr 


( 2 . 12 ) 


where x G C The 7L^(f2)-stability of the projection operator i?™ on the space Xh is immediate 

llVR)rxlU^(n) < ||Vx|U2(n) Vy e X;,. 

Given the POD basis, one can exploit it for model reduction in several different ways. One natural 
choice is to use a Galerkin approach, which yields the following reduced-order formulation: with = 
R'hVh G AT™, find C/” £ X™, n = 1, 2,..., N such that 


+ (VG” , = (/(t„), G X^, 

or equivalently with Cq = r(2 — a). 


(2.13) 


6o(C/;^,:^™)+C„r“(Vt/;^, V:^,„) = &„_i([/° ,¥>^)+E(^j-i-^j)(C^r'><Pm)+c„T“(/(t„),^™) G xr 

i=i 

The existence and uniqueness of the POD approximation follows directly by an energy argument 

(see Section below). In the Galerkin framework, the stiffness matrix of the reduced-order formulation 
(2.13) is the projection of that of the global one (2.5) into the POD space X™. It is worth mentioning 


that the degree of freedom of the reduced system (2.13) is m, i.e., the number of POD basis functions in 
X^, which is usually much smaller than that of (2.5), i.e., the number of finite element basis functions. 


This shows clearly the enormous gain in the computational complexity and storage requirement of the 
proposed scheme. 


3 Error analysis 


In this part, we provide a complete error analysis of the proposed scheme (2.13). The discretization 


error consists of three sources: the spatial discretization, temporal discretization and POD approximation. 
It is known that the semidiscrete solution satisfies the following nearly optimal error estimate min], 
where the operator A is the negative Laplacian operator —A with a zero Dirichlet boundary condition. 
The log factor in the error estimate is due to the limited smoothing property of the solution operator 
for subdiffusion, and the prefactor for t —)■ 0, reflects the corresponding solution singularity. 


Theorem 3.1. Let u he the solution of problem (I.l) with £ Lf{Ll), 0 < ct < 1, and f £ 

L°°(0, T; L^(D)), and Uh he the solution of problem (2.2) with Vh = PhV and ft = Phf- Then there 
holds with ih = I log h\ 

\\uit) - Uhit)\\L^(n) < ch'^llh A\l^(0.) + ||/llL“(0,T;L2(n)) 


6 
















Below we derive the errors due to the temporal approximation and the POD approximation that are 
expressed in terms of the data regularity directly. The main novel ingredient in the convergence analysis 
is to establish a suitable stability result for the LI time stepping under realistic assumptions on the data 
regularity. To this end, we shall develop a novel energy argument, based on the monotonicity of a suitable 
quadrature rule. 


3.1 Error analysis of the LI scheme 

Now we develop a novel energy argument for analyzing the LI approximation. We begin with a 
weighted inequality for the weights {bj}, which is crucial for establishing the monotonicity of the quadra¬ 
ture below. 


Lemma 3.1. Let {bj} be defined by (2.4). Then for j = 2,... ,n — 1, there holds 

{j - + {n- < {n + . 

Proof. Using the definition of the weights bj, the assertion is equivalent to: for all j = 2,..., n — 1: 
0' - 1 + - 1) - (n (1 -f -n + j'j {j + t)~°‘dt < 0, 


that is, 


9{t) 


(j - i + t^if + ty 


-dt < 0, 


where the function g : [0,1] —)■ M is defined by g{t) 
with its g'(t) given by 


g'{t) = a 


J - 1 


n(l -I- — n -I- 7 


n + j), 


For a S (0,1), there holds n(l-I-n ^)“ ^ — n-I-j > n^(n-|-1) ^ — n + j = j — n{n + 1) — Hence 

we deduce g'{t) < 0 on the interval [0,1]. It suffices to show that g{0) < 0. Obviously, 


5(0) = {j - inU - - J + n (1 - (1 + n-i)“-i)). 

I 


The term I in the bracket can be rewritten as 

I = J ((1 - - 1) + n (1 - (1 + . 

We claim that the function g{j) = _) (l — (1 — is monotonically decreasing in j. To see this, 

let hft) : (0,1) —>■ K, with h{t) = t~^{l — (1 — tf~°‘). Then h'{t) = —t~^ (1 — (1 — — at)). Next 

consider the function h(t) : (0,1) —>■ K, with h(t) = (1 — t)“. Then hft) = —a(l — and h''{t) = 

(a —l)a(l —< 0, namely, the function h is concave. Then the concavity implies h{t) < h{0) + h'{0)t, 
which gives (1 — t)“ < 1 — at. Consequently, h'{t) > —t“^(l — (1 — at)~^{l — at)) > 0, and hence h 
is monotonically increasing, and the monotonicity of the function g{j) follows. Hence, by the trivial 
inequality (n — \)/n < n/{n + 1), we have 

I < n((l - - 1) + n(l - (1 + 

= n((l - - (1 - (n + < 0, 

which concludes the proof of the lemma. □ 

Now we give an important monotonicity relation of a weighted rectangular quadrature approximation. 
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Theorem 3.2. Let the function / : [0,1] —>■ K he convex and nonnegative with /(O) = 0, and a S (0,1). 
For any n S N, let Xj = ^, j = 0,... ,n, and yj = j = 0,...,n + 1. Then there holds 


n— 1 n 

bjf{xj) < (n + 1)“-^ Y hfiVi)- 
i=o 0=0 


Proof. First we observe the trivial inequalities < ^+T’ Vj ^ ^ l/i+ij for j = 1,..., n — 1. 

There also holds the trivial identity 


j n-j j i j + l 

Xj := - = -- H-- 

n n n + 1 nn +1 

Now by the convexity of the function /, we deduce 

n-j 


n — J J 

- -Vo + -Vo+i- 

n n 


fiXj) = f 


n-J 1 , . 

- Vo + -%+i < 

n n ’ 


/(%) +-/(2/j+i), j = 


With the assumption /(O) = 0, it suffices to consider j > 1 in the sum. Hence 


n 


OL— 1 




i=i 


i=i 




To show the desired assertion, we consider the following three cases separately, first, last and middle 
terms. For the first term, in view of the nonnegativity of the function /, it suffices to show < 

(n + which however follows from a G (0,1) and 


(n + 1) 


i-a a-i^i-l /n + l\ n-1 / - 1 


For the last term, we have 


n — 1 


i_„,n-l 1 /n-1 


6„_i = n“-i(ni-“ - (n - l)i-“)-= 1- 


n \ n 


< 1 . 


2-0! 


and meanwhile 


(n + l)“-^6„ = (n + l)“-^((n + 1)^-“ - n^-") = 1 - 


n + 1 


1-0 


Hence, it suffices to show n~^ + (1 — n“^)^““ — (1 — (n + 1)“^)^““ > 0 for n > 1. Let g : [0,1] —)■ M by 
g{t) = t + (1 — t)^““ — (1 + Then g{0) = 0, and g'(t) = 1 — (2 — a)(l — t)^““ — (a — 1)(1 + 

Clearly g'{0) = 0 and further g"{t) = (2 — a)(l — a)((l — t)~°‘ — (1 + t)°‘~^) > 0, which in particular 
implies g'{t) > 0 on the interval [0,1]. To conclude the proof, it suffices to show the inequality for the 
middle terms, i.e., for j = 2,..., n — 1 

n“-i^^6,_i +n“-i'^6, < (n + 1)"-^^,, 

n n 

which however is already shown in Lemma |3.1[ □ 

The following result is a direct corollary from Theorem |3.2[ and it will play a crucial role in establishing 
the stability result in Theorem |3.3| below. 























Lemma 3.2. For any a S (0,1), let bj be defined in (2.4|. Then for any n S N, there holds 


- 6,)(n + 1 - jT-i < (n + 
i=i 

Proof. Consider the function /(x) = (1 — x)°‘~^ — 1. Then it satisfies /(x) > 0, /(O) = 0, and also 
/"(x) > 0, i.e., convex. Hence, by Theorem 3.2 we have 

n—1 n 

^ b,{{l - jn-^r-^ - 1) < (n + l)“-i ^ b, ((1 - j[n + - l). (3.1) 

3=0 3=1 

Meanwhile, it can be verified directly that for all n G N+j 1 bj = (1-a) fg X “ dx = 1, i.e., 

^a-i hj = (n + 1)““^ 1^7=0 bj- Plugging the preceding identity into (3.1) yields 


n—1 


—1 


^'^bj{l-jn i)“ <(n + l)“ ^'^bj {1- j{n+l) i) 

3=0 j=0 


a— 1 


which upon rearranging terms gives the desired assertion. 


□ 


Next we give an important L^(n) stability result. The stability estimate puts more weights on the 
source term Fj^ as the index k gets close to the current time step n, in a manner analogous to the 
continuous problem. 


Theorem 3.3. Let Uf, n = 1,2,..., N, be the solution of the fully discrete scheme (2.5). Then with 
Ca = r(2 — a), for n = 1,2,..., N, we have the following stability estimate 


\\Uf\\L'^{n) < \\vh\\L^{n) + CaT°‘ - fc)“ 


(3.2) 


fe =0 


Proof. We show the assertion by mathematical induction. First we consider the case n = 1. Multiplying 
both sides of (2.5) by Ul and integrating over the domain ft yield 

= {UlUl) + c^r^{Fl,Ul). 

Then the Cauchy-Schwartz inequality and Young’s inequality give 

\\Uh\\L^{n) < l|C^°llL=(n) + Ca'r“||F’^||i2(Q). 

Now assume the estimate holds up to some n > 1. A similar argument yields 

n 

||C/r'||L=(0) < 6„||C/^°|k2(n) +^(6,-1 -6,)||t/"+i-^'|U2(n) +c„r“||F-+i|U2(n) 

3 = 1 

n , 


3 = 1 




-fe + l| 




= \\u? 


hllLHQ) + CaT^^^ibj-i - bj)^{n + 1 - j - fc)“ -hCaT“||F"+i||i2(n). 

j—1 k—0 
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Then by changing the order of summation and applying Lemma |3.2| we have 

n n—j n—1 n—k 


- 6,)(n+1 - j - kr-^ 




/c=0 


k^O 

n—1 




<Y.\\PH^'hHn){n + l-kr 


fe=0 


and consequently 


n—1 

||C/r'||L=(a) < iKhnn) + c.r“ ^ " ^)“-'ll^’r'llL^(0) + c„r“||F"+i|U.(a) 

n 

= l|C/.°l|L=(0)+c„T“^(n + l-fcr-i||Fyi|U.(a), 

fc=0 


which completes the induction step and the desired assertion follows. □ 

The next lemma gives one useful estimate for bounding the local truncation error. 

Lemma 3.3. For any 5 S (0, a], there exists a constant c > 0, independent of n, such that for all n>2 

n—1 

((« - -{n-k- 1)1-“) k^-^ < c{n - l)-“. 

k=l 

Proof. The case n = 2 is trivial, and we consider only n > 3. Let dk = ((n — /c)i-“ — {n—k— l)i-“)fc^-^. 
First, we observe that for fc = 1 


di = (n—1)1 “ —(n —2)1 “ = (1 — a) J {n — s) ds 
<c(n — 2)-“ < c((n — l)/3)-“ < c(n — 1)““. 


The sum of the remaining terms can be bounded directly by 

n-l n-l .k+1 n-1 


dk = y](l ~ a)k^~'^ / (n — s)-“ ds < / (n — s)~^{s — l)"^”^ ds 

k—2 k—2 ^ k—2 ^ 


rn—1 


= cy (n — s) “(s —1)° ds = c J (n — s — 1) “s” ^ ds 

n-l 

\-a„5-2j„ I „ / /„ „ i\-a„i5-2. 


/ n — 1 

(n — s — l)-“s‘'-^ ds := I + II. 


Then the desired result follows from 


and 


l<cj (n —s —1) “s° ^ ds < c(n — 1) “y s° ^ds<c(n—1)' 


/ n — 1 

(n - s - l)-“ ds < c(n - l)^-“-i < c(n - 1)"“. 


□ 
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Next we derive an error bound on the local truncation error r„ defined by 

rn = \\dtUh{tn) - d'^Uh{tn)\\L^(n), n = 1,2, (3.3) 

In view of Theorems A.l and |A.2| in the appendix, we make the following temporal regularity assumption. 
Assumption 3.1. The solution u satisfies the following smoothing properties 

||M( 0 llL 2 (n) < c and ||5™M(t)||L2(n) < 
where i5 > 0 and the integer m > 1. 


Remark 3.1. By Theorems \A . 1\ and \A.^ the regularity condition in Assumption \3.1\ holds with S = aa, 
a G (0,1], for initial data v G D{A‘^) and source term f G lT^’°°(0,T;L^(n)). Under these conditions, 
Assumption \3. 1\ holds also for the semidiscrete Galerkin approximation Uh, with a constant c independent 
ofh. 


Lemma 3.4. Let Assumption 3.1 hold, and rn he the local truncation error defined by (3.3). Then 

-5—cx. 


Tr,. < 


CT" “ if n = 1, 

c{n — ifn>2. 


Proof. Using Assumption 3.1 for n = 1, we have the following estimate (with c'^ = l/r(2 — a)) 

/ iT-s)-°‘ [ (u'y,{s) - u'h{y)) dy ds 
Jo Jo 


ri < c'„T ^ 


L2(n) 


<C'„T ^ 




f (t-s) °‘ [ ||<(s)||i 2 (o) + ||<(?/)||i 2 (n)d?/ds 
Jo Jo 

Jo Jo 


(3.4) 


Now we consider the case n > 2. Then 


n—1 


/ II \-a( If \ ’’^hitk+l) — Uh{k) 

= CallZ^/ [tn-s) M^(s)--- 


fc =0 


ds||L2(a) 


<c^l!/ (t„ - s)"“( <(s)- 

fc=0 ^ 


Uhftk+i) - Uh{k) 


) n—1 

k^O 


The first term can be bounded using Assumption 3.1 and the argument for (3.4) as 


rn,o<c (t„-s) “||<(s)||L 2 (f 2 )ds + cr / (t„ - s) 


\uh{y)\\L^n) dyds 


(3.5) 


<c{tn — ti) “ / s^ ^ ds + CT^ ^ I {tn — s) ds < c{n — 1) 

Jo Jo 

Next we derive estimates for rn,k^ A: = 1,2, ...,n — 1. To this end, we use the identity 

Uh { tk + i ) - Uh { k ) 1 /■‘'“+1 .. 1 /■‘'■+1 r 


u'his) - ■ 


r^k+1 ^ ptk+i ps 

/ u'i^is)-u'f,{y)dy=- / ul{z)dzdy 

Jtk ^ Jtk Jy 


and apply Assumption 3.1 such that \\u'/,^{z)\\]^ 2 (^q^ < ct^f, ^ with c independent of t and h to deduce 


4(5) - 


Uh{tk+i) - Uh{k) 


< 


L2(n) 


rtk+l |•ma■x{s,y) 
Itk v/min(s,y) 


K(z)IIl2(0) dzdy < cTt^ 


6-2 
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Thus we obtain 


'^n.k ^ 


cTtl~‘^ / - s)~“ ds = ((n - - {n-k- 1)^““) 

Jtk 

= cT^-°‘k^-'^ {{n - kf-^ -{n-k- l)i-“) . 


Then by Lemma 3.3 we deduce 


rn,k < CT^~°' Y ((^ “ -{n-k- 1)^-“) < CT^-°‘{n - l)’ 


fe=i 


k=l 


This together with (3.51 yields the desired estimate and hence completes the proof. 


□ 


Next we derive the error estimate = Uh{tn) — UJ^, n = 1,2,..., N. First, we observe that the nodal 
error eJJ satisfies = 0 and the following error equation 

d?el + Ahcl = d^Uh{tn) - d^Uh{tn)- 

The next theorem gives an optimal (uniform in time t) error estimate for the fully discrete scheme 


(1^. 

Theorem 3.4. Assume f G W^'°°{0, T; {12)) and v G D{A‘^), with 0 < cr < 1. Let Uh and be the 
solutions of problems (12.21) and (2.5), respectively. Then there holds 


\\Uh{tn) - U]f\\L^{n) < (|!^'^^^||L 2 (n) + ||/||lV2,c=o(o,T;L2(n))) • 

Proof. By Theorem |3.3| and Lemma |3.4[ with 6 = aa, we have 

n—1 


Wh{tn) — b^^||L 2 (o) < cr“ ^(n — k)^ ^\\d^Uh{tk-\-i) — d^Uh{tk-\-i)\\L^{n) 

k^O 

< (||A^z;||L2(n) + ||/||iy2,oo(0,T;L2(n))) 

^ k^l ^ 


Then the following uniform bound 

n— 1 


” ^ 1 ’t ^ /h\~°‘ 

Yin - k)<^-^k-- = n g (' - n) (n) ^ ^ 


yields the desired estimate. 


□ 


Last, we can state an error estimate on the fully discrete approximation Uf, which follows from 
Theorems |3.1|and|3.4|by the triangle inequality. 


Theorem 3.5. Assume f G W'^’°°{0,T-, L^{12)) and v G D{A°'), with 0 < tr < 1. Let u and Uf be the 
solutions of problems (1.1) and (2.5), respectively. Then with ih = |log/i|, there holds 

\\u{tn) - Uf;\\L 2 (n) < + c{h'^£l + t'^“)||/|| w2,oo(o_T;L2(n))- 


3.2 Error analysis of the POD approximation 

Next we derive the error estimates for the POD approximation U)f. First we recall an approximation 
property of the Ritz projection operator R™ defined in (2.12) within the ensemble m Lemma 3 and 
Corrolary 3]. 
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Lemma 3.5. For every m = 1, ...,r, the Ritz projection operator i?™ satisfies 


1 ^ 


and 


1 

N 


n—1 


N 


j—rn-\-l 


A, 


5] (||V(t/," - i?rc/.")lli^(a) + l|V(a“C7,- - 9“i?rc/.”)Hi.(n)) < A, 

n—1 


where {Aj}J^]^ and {^j}j^i denote the eigenvalues of K and K defined in (2.8| and (2.10), respectively. 

Now we can give the error estimate for the POD approximation Uff for smooth problem data. The 
result indicates that the error incurred by using the POD basis in place of the full Galerkin FEM basis is 
determined by the eigenvalues corresponding to the eigenfunctions that are not included in constructing 
the POD approximation. In particular, if the eigenvalues of the correlation matrix decay rapidly, then a 
small number of POD basis functions in the Galerkin POD scheme (2.13) suffice the desired accuracy. 


Theorem 3.6. Let u and Uff be the solutions of (1.1) and (2.13), respectively, and suppose that v G 
D[A), and f G T; L^(0)). Then there holds 


and 


^ ^ ll^(^n) — G^|li2(n) < ct f ^ A^j 

n—1 ^ j—m-\-l ' 

^ ^ ~ < Ct ("t^“ + + h ^ ^ Aj V 

n=l i—m-\-l ' 


(3.6) 


(3.7) 


where {^j}’j=i and denote the eigenvalues of K and K defined in (2.8) and (2.10), respectively. 

Proof. We split the error = u{tn) — Uff into 

e” = Htn) - UJ:) + (UJ: - Uff ), 
and the first term can be bounded using Theorem |3.51 i.e., 

N 

N 


1 « 

■/V 51 ~ ^hWl^in) ^ c (r^“ + 


Hence it suffices to establish a bound for the second term UJ^ — Uff. Now we consider the splitting 

u;: -u:f = {ufi - Rffu^) + {rtuj: - ud ■.= p-+ r. 

Then Lemma [3.5] yields the following bound on p" as 

N r 


1 ^ 

N ^ 


P"lli2(n) ^ c ^ Xj 

j=m-\-l 


and 


-E 

N ^ 


P"lli 2 (o) < ch ^ 


E 

j=m-\-l 


(3.8) 


for the and L^(fl)-POD basis, respectively. Next we derive an estimate on the component 0". 

Using (2.13), the definition of the Ritz projection operator i?™, and the fact that ipm G XJfi C Xh, we 
have 

(a“r, + (vr, = {d^KUff, g^m) + vp™) - ( 5 “c/;^, p™) - (vc/”, v^p^) 

= {d:fR]:UY g^m) + mj:, Vipm) - {f{tn), Tm) 

= - UY), Tm) = -(a“p”, g,m) 
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and 0*^ = 0. The stability result in Theorem 3.3 yields 


n—1 


|| 0 ”||L^(n) < cr“ J 2 in- 


k^O 


Appealing to Young’s inequality for the Laplace type discrete convolution [Tj Theorem 20.18], i.e., 


N 


N 


N 




(3.9) 


n—O \fc=0 


\n—0 


n—0 


we deduce 

N / 71—1 


N yn-1 ^2 y N K 2 N N 

n—1 ^ fc=0 ^ ^ n—1 ^ 71—1 71—1 


fc=0 

Then by Lemma 3.5 we have 

N 

N 


— y I 

N ^ ' 

n—1 


n\\2 

\\L^m ^ 


< Y WdrP^Wlnn) = ^ Y WdrP'^Wl^n) < ct 

n—1 n—1 j—7n-\-l 


Y A. 


Likewise, for the L^(n)-POD basis, we deduce 
1 ^ 


— y ||0 

N ^ " 

n—1 

This completes the proof of the theorem. 


n—1 n—1 


<CTh'^ Y Aj- 

j—rn-\-l 


□ 


The error estimate in Theorem 3.6 covers only smooth initial data v S D[A). In the case of nonsmooth 
initial data v G D{A'^), 0 < a < 1, one can derive an analogous error estimate; see the following remark. 

We note that the regularity of problem data (or solution) does not enter the error estimate due to the 
POD approximation directly. Hence, in principle, the approach is capable of handling nonsmooth problem 
data, if the solution singularity is built-in in the ensemble of snapshots and thus captured by the POD 
basis directly. 

Remark 3.2. We comment on nonsmooth problem data. Consider the Hl{^) POD for f G 1T^’“(0, T; L^(r2)) 
and nonsmooth initial data v G D{A'^), 0 < cr < 1. Then in view of Theorem S.5, we have 

N N 

- y iiu(t„) - + h%- y 


n—1 


n—1 


Meanwhile, the summation can be bounded as 

1 ^ _-2q(1-ct) N _-2a(l-cr) rN 

1 y- ,-2a(l-a) _ ^ 2a(l-^)^^ ' 

TV ^ " N 2^ - 


n—1 n—1 

where the constant ia,a,T is given by 


N 


S '^^ds < CT^a.cr,! 



), a(l-(T)>l/2, 


a(l - a) = 1/2, 
a(l-CT) < 1/2. 
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Consequently, by repeating the arguments in Theorem \3.6\ we obtain the following error estimate for the 
POD approximation {U^} (with the HqITI) POD basis) 


n—1 ^ ^ 

and a similar error estimate holds for the L^(r2) POD basis. Interestingly, for the case a{l — a) < 1/2, 
the error estimate in the space remains uniform with respect to the time step size r. 

Last we briefly comment on the case when the FDQs are not included in the snapshots. 

Remark 3.3. In our construction of the POD basis, we have included the FDQs in the snapshots. When 


the FDQs dfUf), n = 1,2,..., N, are not contained in the snapshot set, the error formula (2.9) for 


POD basis becomes 


1 


N 


+ 1 




n—0 


i=i 




Further for the FDQs we have 

N 

N 


.. N m m 

- E WdrUj) - E(v9“c/i/, =- E 


\<p) 


j = l n=l j=l 

Let = U)) — ,'V'ipj)'f’j. By the monotonicity of the weights {bj}, we have 

n—1 

— _ n ^ ~ / />o n 

drUn 


< CaT boWU+ bn-l\\Ui^\\H^(n) + E(^l-1 “ ^j)\\Uh 

, n S 2 

<CaT + CaT ’ 


J=0 


with gj = bj-i — bj and 6_i = 2. Then by Young’s inequality for discrete convolution, cf. (3.9), we arrive 
at 


.. N , n •.2../N..2N N 

n—1 ^ j—0 ' ^ n—0 '' n—0 n—0 




Meanwhile, by the Cauchy-Schwarz inequality, we have 


^ ^ / rn+l \2 I-N+I I ’ */< 1/2, 

E = (1 ~ E ( / t clog TV, if a = 1/2, 

^ * c, if a > 1/2. 


Consequently, there holds 


1 

TV 


N 


n=l j=l 

where the constant la_,T is given by 



< CT{la 

,r\\U, 

1 

if a < 

1/2, 

logy 

if a = 

1/2, 

.-2a+l 

if a > 

1/2. 
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j—rn-\-l 










For a < 1/2, the term involving the initial data [7° is of higher order in comparison with the last term. 
Hence the error for HQ{fl) Galerkin POD (2.13) (without FDQs in the snapshots) can be bounded by 


N 


n—1 


- ^ ||u(t„) - < CT + hX+io.AU°H - 


+ T 


-2a 


i=i 


E A,). 

j=m+l 


In comparison with the error estimate (3.6) with FDQs from Theorem \ 3.6\ this estimate contains an 
extra factor and an approximation error of the initial data Vh (within the POD basis XJf). For the 
fractional order a —^ 1, the factor recovers that for the classical diffusion equation im- 


4 Numerical results 

Now we present numerical results to verify the convergence theory in Section and the efficiency of 
the proposed Galerkin-Ll-POD scheme. 


4.1 Numerical results for one-dimensional examples 

First we present numerical results for one-dimensional examples to verify the convergence analysis in 
Section We consider the subdiffusion model in the following two cases: 

(a) n = (0,1), v = xil-x)G D{A), and fix,t) = e 

(b) VL= (0,l),n = X(o,i/2)(a;) € for e G ( 0 , 1 / 4 ), and/(x,t) = g W^'°°{Q,T-,L^{FL))- 


In the computations, we divide the unit interval VL into M equally spaced subintervals with a mesh size 
h = 1/M. Likewise, we fix the time step size t at r = T/N. 

First we examine the temporal convergence by setting T = 0.1 (the spatial convergence was already 
examined in [lain). We take a small mesh size h = 10 so that the spatial discretization error is 
negligible. The exact solution can be expressed in terms of the Mittag-Leffler function Ea^p{z), cf. (A.l), 
which can be evaluated efficiently by an algorithm developed in [32) . The numerical results by the fully 
discrete scheme (2.5) are given in Table In the table, rate refers to the empirical rate when the time 


step size r halves, and the numbers in the bracket denote the theoretical predictions from Theorem 3.5 
For cases (a) and (b), the empirical rate is and 0(r“/"^), respectively, which agree well with the 

theoretical ones. The convergence rate of the LI scheme improves with the smoothness of the initial data 
V (while keeping the smooth right hand side / fixed) and the increase of the fractional order a, since the 
solution regularity improves accordingly. 


Table 1: The maximum error emax = inaxi<„<Ar ||t7^ — u(tn)\\L^{n) for initial data (a) and (b) with 
T = 0.1, h = 10-3, T = T/N. 


a 

N 

1000 

2000 

4000 

8000 

16000 

32000 

rate 

0.35 

(a) 

2.67e-3 

2.27e-3 

1.90e-3 

1.58e-3 

1.29e-3 

1.05e-3 

« 0.29 (0.35) 


(b) 

2.48e-2 

2.41e-2 

2.29e-2 

2.15e-2 

1.99e-2 

1.82e-2 

« 0.10 (0.09) 

0.5 

(a) 

9.26e-4 

6.73e-4 

4.86e-4 

3.50e-4 

2.51e-4 

1.80e-4 

« 0.48 (0.50) 


(b) 

2.03e-2 

1.81e-2 

1.64e-2 

1.50e-2 

1.37e-2 

1.26e-2 

« 0.13 (0.13) 

0.75 

(a) 

1.82e-4 

1.09e-4 

6.43e-5 

3.77e-5 

2.17e-5 

1.25e-5 

« 0.76 (0.75) 


(b) 

2.52e-2 

2.20e-2 

1.91e-2 

1.64e-2 

1.39e-2 

1.15e-2 

« 0.21 (0.19) 


Next we illustrate the proposed Galerkin-Ll-POD scheme, and the numerical results are given in 
Table for the choice T = 1 and N = 200. Here the average error e and the POD approximation error 
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e™ are defined by 


and 

n—1 n—1 

respectively. Like before, we use the notation ~ and ^ over e™ to denote iLQ(r2)-and L^(fl)-POD basis, 
respectively, and the subscript w to indicate that the snapshots do not contain FDQs. For example, e™ 
and denote the error between the full Galerkin solution UJ^ and the solution of the Galerkin POD 
formulation with m Lfo(O) POD basis functions, with and without FDQs, respectively. For both cases 
(a) and (b), with three or four POD basis functions, the POD approximation error falls below the error 
due to temporal discretization, and the convergence is relatively independent of the fractional order a. 
The fast convergence of the Galerkin POD scheme is also expected from the exponential decay of the 
eigenvalues of the correlation matrix, cf. Fig. [l] Further, the inclusion of FDQs does not affect much the 
POD approximation error, with their errors within a factor of ten, even though their presence improves 
the apparent theoretical convergence rates, cf. Theorem |3.6| and Remark |3.3| The effect seems to be 
compensated by the smaller eigenvalues, cf. Fig. These observations show the efficiency of the Galerkin 
POD scheme, which has only a degree of freedom of three or four at each time level, compared with one 
thousand for the standard Galerkin FEM. 

For case (b), the Galerkin POD scheme requires slightly more POD basis functions in order to reach 
the same level of the accuracy. This is expected, since for nonsmooth data v, it can only be accurately 
described by more Fourier modes, and all these modes persist in the dynamics due to the “slow” decay of 
subdiffusion. Hence the solution manifold may exhibit richer structure than case (a), and consequently, 
more POD basis functions are needed to accurately capture the dynamics. However, the eigenvalues in 
the nonsmooth case decays also exponentially, cf. Fig. Hence, the proposed scheme also works well 
with low regularity data. 

The efhciency of the proposed scheme relies crucially on constructing “good” POD basis. To this end, 
we present the first five POD basis functions for case (b) in Fig. The iLQ(O)- and LF’iTl) POD basis 
take very different shapes: for the Hl{Vl) POD, the first basis function captures the singularity (caused 
by the discontinuous initial data), whereas the higher POD modes are very smooth. In contrast, for 
the L^(0) POD, all the first five POD basis functions contain singularities (in the middle of the interval 
as well as oscillations around the end points). Namely, the i?o(0) POD seems to better aggregate the 
solution singularity (actually into one single POD basis). Nonetheless, the L^(0) and i?o(0) POD-basis 
exhibit quite similar approximation property, and thus can provide equally good approximations of the 
solution manifold, cf. Table 

4.2 Numerical results for one two-dimensional example 

Now we present numerical results for the following two-dimensional example: 

(c) O = (-1,1)^\([0, l]x[-l,0]), v{xx,X 2 ) = a:i(l-l-a;i)(l-a:i) sin(27ra;2), f{xi,X 2 ,t) = ™(7i-a:2) g 

and r= 1. 

In the computations, we divide the L-shaped domain Q into a triangulation with a degree of freedom 
lO'^, and fix the time step size r at t = T/200. A reentrant corner with an angle ui G (tt, 27r) induces a 
singularity associated with the corresponding stationary Poisson’s problem [5] . In example (c), the angle 
ui = 37r/2, and the reentrant corner gives rise to a singularity near the origin with a leading term of the 
form sin(26*/3) in polar coordinates. Hence, we refine the mesh adaptively using the bisection rule 
[IZl Section 4.1]. We compute the reference solution on a more refined mesh with 2 x 10^ and r = 1/1000. 

The numerical results are shown in Table The POD scheme exhibits a fast convergence, and the 
error decreases steadily with the increase of the number m of POD basis functions. In particular, five or 
six POD basis functions suffice to resolve the solution manifold to an accuracy 0(10“®), which clearly 
shows the efficiency of the Galerkin POD scheme, when compared with the standard Galerkin FEM. 
The fast convergence follows also from the exponential decay of the eigenvalues of the correlation matrix. 
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Table 2: The numerical results of the Galerkin POD for cases (a) and (b) with T = 1, /i = 10 TV = 200, 
and with m POD basis functions. 


a 

case 

m 

e 


'pn 

gm 



(a) 

3 

1.82e-7 

9.34e-12 

3.03e-12 

9.45e-12 

3.02e-12 



4 

1.82e-7 

4.72e-13 

3.71e-14 

4.83e-13 

3.19e-14 

0.3 

(b) 

3 

3.83e-6 

4.65e-6 

3.59e-6 

4.36e-6 

3.60e-6 



4 

3.83e-6 

2.73e-9 

2.41e-9 

2.73e-9 

2.41e-9 


(a) 

3 

4.46e-7 

l.Ole-10 

6.25e-12 

l.lle-10 

6.22e-12 



4 

4.46e-7 

5.33e-13 

8.87e-14 

5.41e-13 

8.28e-14 

0.5 

(b) 

3 

1.70e-5 

1.81e-5 

6.70e-6 

1.59e-5 

7.08e-6 



4 

1.70e-5 

3.67e-8 

6.70e-9 

3.43e-8 

6.69e-9 


(a) 

3 

2.89e-7 

4.70e-10 

1.35e-ll 

4.98e-10 

1.34e-ll 



4 

2.89e-7 

1.33e-12 

1.85e-13 

1.29e-12 

1.81e-13 

0.7 

(b) 

4 

2.80e-5 

2.51e-5 

1.83e-7 

1.45e-5 

1.78e-7 



5 

2.80e-5 

2.49e-8 

5.00e-9 

2.42e-8 

4.99e-9 



(a) case (a), a = 0.3 (b) case (a), a = 0.5 (c) case (a), a = 0.7 





(d) case (b), a = 0.3 


(e) case (b), a = 0.5 


(f) case (b), a = 0.7 


Figure 1: The decay of eigenvalues of the correlation matrix in the ID problem with a = 0.3, 0.5 and 0.7. 
Here, A„, A™, A„, and denote eigenvalues of correlation matrix for i?o(D) POD basis with or without 
FDQs and L^(D) POD basis with or without FDQs, respectively. 


cf. Fig. The decay rate of the spectrum is almost identical for the L^(D) and POD basis, 

and independent of the presence of the FDQs. Hence, the presence of geometrical singularities in the 
domain does not influence the efficiency of the Galerkin-L 1-POD scheme. Interestingly, we observe that 
with the increase of the fractional order a, the error increases slightly, which awaits further theoretical 
justification. 
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Figure 2: The first five POD basis functions, in the and L^(0) norms for case (b), with FDQs 

included in the basis construction. 


Table 3: The numerical results of the Galerkin POD for case (c), with T = 1, N = 200 and with m POD 
basis functions. 


a 

m 

e 

~m 

-pn 

e™ 

Cm 

0.3 

5 

7.67e-7 

5.36e-10 

3.33e-10 

5.17e-10 

3.32e-10 


6 

7.67e-7 

6.40e-12 

5.49e-12 

6.39e-12 

5.48e-12 

0.5 

5 

4.75e-6 

2.08e-8 

8.23e-9 

1.96e-8 

8.18e-9 


6 

4.75e-6 

1.62e-10 

4.82e-ll 

1.44e-10 

4.79e-ll 

0.7 

6 

l.Ole-5 

2.05e-8 

1.36e-9 

1.38e-8 

1.27e-9 


7 

l.Ole-5 

9.11e-10 

1.17e-10 

6.09e-10 

l.lle-10 





Figure 3: The decay of the eigenvalues of the correlation matrix for case (c) (2D problem on an L-shaped 
domain), with a = 0.3, 0.5 and 0.7. Here A„, A“, A„, and AJC denote eigenvalues of correlation matrix 
for 7 Jq(H) pod basis with or without FDQs and L^(D) POD basis with or without FDQs, respectively. 
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4.3 Numerical results for a perturbed problem 

Last, we illustrate the proposed Galerkin POD scheme with a perturbed problem, where the snapshots 
are generated using a problem setting different from the one of interest, as typically occurs in optimal 
control and inverse problems. Let Sn{x) = n{2 cosh^{nx))~^ be an approximate Dirac delta function. 

(d) On the domain D is O = (0,1)^, we consider the following problem: 

d“u — Am + qu = f in D 

with q{xi,X2) = 1 + cos(7ra:i) sin(27ra:2), f{xi,X2,t) = 62{xi — \)52{x2 — and v{xi,X2) = 

a:i(l — Si) sin(27ra:2) and T = 1. However, the snapshots are generated using a perturbed source 
term f{xi,X 2 ,t) = (5io(a:i - l)Sioix 2 - \). 

In our computation, we divide the sides of the domain D into 100 equal subintervals, each of length 
10“^, thus dividing H into lO'^ small squares, and obtain a uniform triangulation by connecting parallel 
diagonals of each small square. The time step size r is fixed as t — T /200. 


Table 4: The numerical results of the Galerkin POD for case (d), with T = 1, N = 200 and with m POD 
basis functions. 


a 

m 



e™ 


0.3 

4 

4.63e-7 

4.64e-7 

4.63e-7 

4.64e-7 


5 

3.32e-7 

4.50e-7 

3.2Ie-7 

3.34e-7 

0.5 

4 

4.47e-7 

4.52e-7 

4.47e-7 

4.53e-7 


5 

3.50e-7 

3.46e-7 

3.50e-8 

3.45e-7 

0.4 

4 

4.I2e-7 

4.32e-7 

4.I2e-7 

4.32e-7 


5 

3.81e-7 

3.7Ie-7 

3.80e-7 

3.71e-7 


Since the snapshots are generated from a perturbed problem, the error estimates in Theorem |3.6| do 
not apply directly. Nonetheless, one can still observe a fast decay of the POD approximation error, and 
with four to five POD basis functions, the error is already much smaller than the LI time stepping, 
cf. Table for both L^(r2)- and iLQ(r2)-POD basis and with/without FDQs. The high efficiency of 
the proposed scheme is attributed to the intrinsic low-dimensionality of the solution manifold, which is 
fully captured by the snapshots generated from the perturbed problem. This is also expected from the 
fast decay of the eigenvalues of the correlation matrix (from the perturbed problem) in Fig. The 
solution profiles and corresponding errors are shown in Fig. This example shows clearly the potential 
of the proposed approach for solving related inverse problems and optimal control, where many analogous 
forward problems have to be solved. 

5 Concluding remarks 

In this work, we have developed an efficient Galerkin-Ll-POD scheme for solving the subdiffusion 
problem, by coupling the Galerkin finite element method, LI time stepping and proper orthogonal de¬ 
composition. It realizes the computational efficiency by constructing an effective reduced-order model 
using POD, often with a very small degree of freedom. We provided a complete error analysis of the 
scheme, and derived optimal error estimates due to spatial discretization, temporal discretization and 
POD approximation. This is achieved by developing a novel energy argument for LI time stepping. 
The extensive numerical experiments fully confirmed the convergence analysis and the efficiency and 
robustness of the scheme. 

The work represents only a first step towards effective model reduction strategies for fractional differ¬ 
ential equations. The choice of the three components in the proposed scheme is not unique. Alternatively, 
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Figure 4: The decay of the eigenvalues of the correlation matrix for case (d) with a = 0.3, 0.5 and 0.7. 
Here A„, AJJ", A„, and A“ denote eigenvalues of correlation matrix for Ho(fi) POD basis with or without 
FDQs and L^(0) POD basis with or without FDQs, respectively. 





(a) exact solution, a = 0.3 


(b) POD solution, a = 0.3 


(c) error, a = 0.3 



(d) exact solution, a = 0.5 


(e) POD solution, a = 0.5 


(f) error, a = 0.5 



(g) exact solution, a = 0.7 


(h) POD solution, a = 0.7 


(i) error, a = 0.7 






Figure 5: Exact and numerical solutions at T = 1 for case (d), where the POD solutions are obtained 
using POD basis with the FDQs. 


one may employ finite difference methods or spectral methods instead of the finite element method, and 
convolution quadrature type schemes instead of the LI time scheme. The overall framework extends 
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straightforwardly to these alternative choices, even though the convergence analysis will differ. Further, 
it is of much interest to extend the proposed scheme to more complex models, e.g., the multi-term model 
and the diffusion-wave model. 
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A Regularity theory for problem (1.1) 


Now we describe temporal regularity results of problem which plays an important role in the 
convergence analysis. Let {(Aj, be the eigenvalue pairs of the negative Laplacian A = —A with a 

homogeneous Dirichlet boundary condition, where the set forms an orthonormal basis in 

Then by the standard separation of variable technique, we deduce that the solution u can be represented 

by 

u(t) = E(t)v + / E(t — s)f{s)ds, 

Jo 

where the solution operators E{t) and E{t) are given by 


E{t)ij} = '^Ea^i{-Xjt°‘){i’,(pj)(Pj and E{t)i{j = ^Ea,a{-Xjt°^){ilj,(pj)(pj, (A.l) 

f=i 

respectively. Here the Mittag-Leffler function Ea^p{z), a > 0, /3 G M, is defined by [TSJ pp. 42] Ea,p{z) = 
r{ka+p) ■ "bbe following relations hold (see [211 Lemma 3.2] and [HI pp. 43, eq. (1.8.28)] for proofs). 

Lemma A.l. Let a G (0,1), and /3 G M. The Mittag-Leffler function Ea^piz) satisfies for m > 1 

dm 

— E„.i(-Af“) = -Af“-™E„,„+i_„(-At“) t > 0, 

at"* 

and the following uniform bound on the negative real axis ]R“ holds 

Ea,l3iz) < c{l\z\)~^ VxGM”. 


Now we can state the temporal regularity for the homogeneous problem. 
Theorem A.l. If v G and f = 0, then 


\\(Itu\\mn) < ct'^^-"^\\A-v\\L2in), 

where if a £ (0,1], m > 1, and if a = 0, 0 < m < 2. 

Proof. The case tr = 0 has been shown m Corollary 2.6]. Eor tr G (0,1], by Lemma [A. 1[ we have 

oo 

i-1 


/\ j.ClL\2 — 2<7 ^ 

j (L -t Ajl ) 

where the last inequality follows from the inequality sup^ (Ajt“)^“^'^/(1 -I- Xjt^^Y < c. 


□ 
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Next we consider the inhomogeneous problem. We shall need the following estimate on E{t) 
Lemma A.2. For any t > Q, we have for x G L^{n) and m > 0 

||ar^(t)xlU.(n)<ct—ixlU=(n). 


Proof. The definition of the operator E in (A.l) and Lemma A.l yield 

OO 

\\drmx\\hin) = J2\t° 


< ct 


oi—m— 1 

h 

2a —2m—2 


£^a,a-m(-Ajt“)| \{x,F])Y 


\{x,Fj)Y = ct 
i=i 


2a —2m —2 ii 


li2(o), 


which completes the proof of the lemma. 

Now we can state the temporal regularity result for the inhomogeneous problem. 
Theorem A.2. If v = 0 and f G (0, T; (12)) with some m G [0,2], then there holds 

ll^r'“IU 2 (n) < '"||/||w™.“(o.T;L 2 (n)), 0 < m < 2. 

Proof. Using the following convolution relation [22l Lemma 5.2] 

tif * g)' = f * g + (tf) *g + f* (tg') 
and Lemma A.2 we deduce that for t G (0, T] 

t™115rulli.(0)<c ^ f \\(t- s)^dlE(t- s)(s‘^f^^\s))\\L.^n)ds 

, ^ -'0 

p+g<m 

p+g<m p+q'Cm 

Since for t G (0,T], we have J2p+q<mt°'^‘^~’^ < CTt°‘~^, the desired assertion follows. 


□ 


(A.3) 


a+q 


□ 


Remark A.l. Theorems \ A.1\ and A.2 show the limited smoothing property of the subdiffusion model 
([T3.- for the homogeneous problem with v G D(A), the first order derivative in time t of the solution u 
exhibits a singularity of the form t°‘~^; and for the inhomogeneous problem with f G W^’°°(0,T; L^(I2)), 
the first-order derivative exhibits a similar singularity, despite the smoothness of f in time. 
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